library("RSvgDevice")
flat <- function (gene.expression.data, data1, data2, indexmethod="intersection", outputplotfile) {

###############  constants used  ################

no.genes <- nrow(data1)

no.clusters.1 <- max(data1[,2])
no.clusters.2 <- max(data2[,2])

clus <- no.clusters.1 + no.clusters.2

no.clusters <- max(no.clusters.1,no.clusters.2)

limit <- max(abs(min(gene.expression.data, na.rm=T)),abs(max(gene.expression.data, na.rm=T)))

time.points <- ncol(gene.expression.data)


################  defining the colours ####################

step <- (floor(100/(ceiling(clus/3))))/100

colours <- matrix(0,1,3)

for (i in 1:ceiling(clus/3))
  {
     for (j in 1:3)

       { 
        try.again <- TRUE
        while (try.again) 
         {
           try.again <- FALSE
           new.colour <- c(ceiling(runif(1,-0.5,1)), ceiling(runif(1,-0.8,1)), ceiling(runif(1,-0.2,1)) ) * step * i 
           for (k in 1:nrow(colours)) {if (identical(colours[k,], new.colour)) {try.again = TRUE}}
         }
         colours <- rbind(colours, new.colour)
       }
   }

colours <- sample(rgb(colours[,1],colours[,2],colours[,3]),3*ceiling(clus/3))


colours.1 <- colours[1:no.clusters.1]
colours.2 <- colours[(no.clusters.1 + 1):clus]


#####################  performing the comparison    ###########################

if (indexmethod == "intersection")

{ comparison <- compare.intersection(data1,data2)}
else
 {
   if ( indexmethod == "simpson") {comparison <- compare.simpson(data1,data2)}
 }

result <- number3 (comparison)

intersection <- compare.intersection(data1,data2)

p <- result$p

###################   grouping the initial clusters in p super-clusters  ###############

new.m2 <- result$m2

visited <- matrix(0,1,no.clusters.2)

index <- 1

for (i in 1:no.clusters.1)
 {
   new <- FALSE

   for (j in 1:no.clusters.2)
      {
        if (visited[1,j] == 0) 
           {
	     if (result$m2[i,j]!=0) 
                {
		  new <- TRUE
                  visited[1,j] <- 1
                  new.m2[,j] <- result$m2[,j] * index
                }
           }  	
      }
   if (new == TRUE) {index <- index + 1}
 }


#####################  relabeling the initial clusters ###########################

new.data1 <- data1
new.data2 <- data2

for (i in 1:no.genes)
    {
      new.data1[i,2] <- max(new.m2[data1[i,2],])
      new.data2[i,2] <- max(new.m2[,data2[i,2]])
    }

############## splitting the screen in p x 6 spots of different width ###############

devSVG(file = outputplotfile, bg = "white", fg = "black", onefile=TRUE, xmlHeader=TRUE)

plot.new()

layout(matrix(c(1:(6 * p)),p,6,byrow=TRUE),width=c(1,2,1,2,1,2))

par(mar=c(1.8,1,1.5,1))

###################################################################


temp <- overlappings(data1,data2,indexmethod)
overlap <- temp$overlapping
Jaccard <- temp$R

for (i in 1:p)  ##### fill in each row

  { 

##### screen 1

    termometer(Jaccard[i,1])
   
##### screen 2

    plot(1:time.points,rep(0,time.points),type="l",main="",ylim=c(-limit,limit),xlab="",ylab="",cex.axis = 0.8, xaxt = "n")
    for (cl in 1:no.clusters.1)
    {
      if (any(new.m2[cl,]== i) ) ## plot this initial cluster
         {
           for (j in 1:no.genes)
            {
              if (data1[j,2] == cl)  {lines(1:time.points,gene.expression.data[j,],col= colours.1[cl])}
            }
          }
    }

##### screen 3

   inter <- matrix(0,1,no.clusters.2)
   for (j in 1:no.clusters.1)
    {
      for (k in 1:no.clusters.2)
         {
           if (new.m2[j,k] == i)  {inter <- inter + intersection[j,]}
         }
    }
    colour.bars (inter,colours.2, width = .5)

##### screen 4

     venn.1(overlap[i,1], overlap[i,2] , overlap[i,3] )

##### screen 5

   inter <- matrix(0,1,no.clusters.1)
   for (k in 1:no.clusters.2)
    {
      for (j in 1:no.clusters.1)
         {
           if (new.m2[j,k] == i)  {inter <- inter + intersection[,k]}
         }
    }
    colour.bars (inter,colours.1, width = .5)


##### screen 6

    plot(1:time.points,rep(0,time.points),type="l",main="",ylim=c(-limit,limit),xlab="",ylab="",cex.axis = 0.8, xaxt = "n")
    for (cl in 1:no.clusters.2)
      {
        if (any(new.m2[,cl]== i) ) ## plot this initial cluster
           {
            for (j in 1:no.genes)
             { 
              if (data2[j,2] == cl)  {lines(1:time.points,gene.expression.data[j,],col=colours.2[cl])}
             }
           }
      }


  }   ### end of the screens

####################################################################

# dev.copy(png, file = "clustering.output.graphics.png")
# dev.off(which = dev.cur())
dev.off()

}

termometer <- function (R, col=2, width=0.25, font.type=2, font.size=1, digits=3) {
  if (R > 1) {stop("the parameter R can't be larger than 1")}

  up <- round(R*10^digits)/10^digits
  down <- 0
  x <- c(0.5-width/2,0.5-width/2,0.5+width/2,0.5+width/2)
  plot.new()
  R.i <- c(down, up, up, down)
                                        # percentage.i <- paste(up, "%")
  frame <- c(0,1,1,0)
  polygon(x,frame, fg="black", col=0, lwd =1)
  polygon(x, R.i, fg="black",col=2, lwd = 1)
                                        # text(0.5, up/2, labels=percentage.i, bg="transparent", col=1, font=font.type, cex=font.size)
  text(0.5, up/2, labels=up, bg="transparent", col=1, font=font.type, cex=font.size)
}

venn.1 <- function (size.A, size.B, intersection, edges = 200, main = NULL, ...) {
#### computing the parameter for the termometer

  union <- size.A + size.B - intersection
  R <- intersection / union

####################################################


  x <- c(size.A, size.B, intersection)
  if (!is.numeric(x) || any(is.na(x) | x <= 0))  {stop("sizes must be positive.")}

  if (intersection > min(size.A, size.B))  {stop("the intersection can't be larger than the sets sizes")}

##### Size of the plot

  pin <- par("pin")
  xlim <- ylim <- c(-1, 1)
  if (pin[1] > pin[2])  {xlim <- (pin[1]/pin[2]) * xlim}
  else {ylim <- (pin[2]/pin[1]) * ylim}
  plot.new()
  plot.window(xlim, ylim, "", asp = 1)

##### Radius proportional to the size of the clusters (=> area) and normalized

  radius.A <- sqrt(size.A/pi)
  radius.B <- sqrt(size.B/pi)

  radius <- c(radius.A, radius.B)

  M <- max(radius)
  radius <- radius/M

  radius.A <- radius[1]
  radius.B <- radius[2]

  area.A <- pi * (radius.A^2)
  area.B <- pi * (radius.B^2)
  area.intersection <- intersection / M^2


##### CENTERS: the the circles are at first drawn centered at (0,0) and (c.B,0)

  epsilon <- 0.01

  continue <- TRUE
  centered <- FALSE

  previous.step <- radius.B
  step <- 0
  previous.direction <- 1

  c.A <- 0
  c.B <- radius.A


  counter <- 0

  if (intersection == min(size.A,size.B)) 
    {
      centered <- TRUE
      continue <- FALSE
      c.A <- 0
      c.B <- 0
    }
  counter <-0 

  while (continue == TRUE)
    {
      a <- (radius.A^2 - radius.B^2 + c.B^2)/(2*c.B)

      y <- (a-c.B)/radius.B
      z <- a/radius.A

      counter <- counter +1

#### to prevent roundoff errors  ####


      if (y > 1) {y <- 1} 
      if (y < -1) {y <- -1}

      if (z > 1) {z <- 1}
      if (z < -1) {z <- -1}


      I.1 <- - (radius.B^2/2) * (acos(y) - pi - (sqrt(1-(cos(2*acos(y)))^2))/2 )
      I.2 <- (radius.A^2/2) * (acos(z) - (sqrt(1-(cos(2*acos(z)))^2))/2)

      area.temp <-  2* (I.1 + I.2)

      if (abs(area.temp - area.intersection) <= epsilon) 
        {
          continue <- FALSE
          direction <- 0
        }
      else
        {

          if (previous.step == radius.B) {step <- previous.step / 2}

          if (area.temp - area.intersection < 0) 
            {
              direction <- -1
              c.B <- c.B - step
            }
          if (area.temp - area.intersection > 0) 
            {
              direction <- 1
              c.B <- c.B + step
            }
        }

      if (direction!=previous.direction) {step <- step/2}
      previous.step<-step
      previous.direction <- direction
    }


  print(c.B)

  print(counter)

  c.B <- c.B/2


################################################################
######################  P I C T U R E  ######################### 
################################################################

  ## Cluster A

  center.A <- (radius.A - radius.B - c.B)/2       
  if (centered) {center.A <- 0}

  t2p <- seq(0, 2*pi, length = edges)
  xc <- center.A+cos(t2p) * radius.A
  yc <- sin(t2p) * radius.A
  polygon(xc, yc, fg="black",col="transparent", lwd = 3)
  lines(x = c(center.A - 1.1 * radius.A, center.A - radius.A), y = c(0,0), lwd =2)
  text(center.A - 1.22 * radius.A, 0, labels=size.A, bg="transparent", col=2, font=2, cex=1)

  ## Cluster B

  center.B <- (c.B + radius.A - radius.B)/2       
  if (centered) {center.B <-0}

  t2p <- seq(0,2*pi, length = edges)
  xc <- center.B + cos(t2p) * radius.B
  yc <- sin(t2p) * radius.B
  polygon(xc, yc, fg="black",col="transparent", lwd = 3)
  lines(x=c(center.B+radius.B, center.B + 1.1 * radius.B), y=c(0,0), lwd = 2)
  text(center.B + 1.22 * radius.B, 0, labels=size.B, bg="transparent", col=2, font=2, cex=1)

  ## Intersection

  center.intersection <- (center.A + radius.A - radius.B + center.B)/2
  if (centered)  {center.intersection <- 0}
  text(center.intersection, 0, labels=intersection, bg="transparent", col=4, font=2, cex=1)


  return(R)

  title(main = main, ...)
  invisible(NULL)
}

compare <- function (datos1,datos2)

{
datos1 <- as.matrix(datos1)
datos2 <- as.matrix(datos2)

num1 <- max(datos1[,2])
num2 <- max(datos2[,2])

tam <- dim(datos1)[1]
int <- matrix(0,num1,num2)


for (i in 1:tam)
{
int[as.numeric(datos1[i,2]),as.numeric(datos2[i,2])] <- int[as.numeric(datos1[i,2]),as.numeric(datos2[i,2])]+1
}


size1 <- matrix(0,1,num1)
for (i in 1:tam)
{size1[as.numeric(datos1[i,2])] <- size1[as.numeric(datos1[i,2])]+1
}

size2 <- matrix(0,1,num2)
for (i in 1:tam)
{size2[as.numeric(datos2[i,2])] <- size2[as.numeric(datos2[i,2])]+1
}

s <- int
for (i in 1:num1)
  {for (j in 1:num2)
    {s[i,j] <- s[i,j]/min(size1[i],size2[j])

    }
  }


return(int, size1,size2,s)
#return(int,s)
#return(int)

}


number3 <- function (matriz) {
################################################
######## MATRIX OF 1'S POINTING OUT THE MAX'S
################################################

  matriz<-as.matrix(matriz)
  r<-dim(matriz)[1]
  s<-dim(matriz)[2]

  m2<-matrix(0,r,s)

  for (i in 1:r)
    {
      maxim <-max(matriz[i,])

      if (maxim >0)
        {
          for (j in 1:s) 
            {
              if (matriz[i,j]==maxim) {m2[i,j]<-1}
            }
        }
    }

  for (j in 1:s)
    {
      maxim<-max(matriz[,j])

      if (maxim>0)
        {
          for (i in 1:r) 
            {
              if (matriz[i,j]==maxim) {m2[i,j]<-1}
            }
        }
    }


#############################################
######### DETERMINING p #####################
#############################################

  m3<-matrix(0,r,1)
  p<-r

  for (i in 1:r)
    {
      for (j in 1:s)
        {
          if (m2[i,j]==1)

            {
              m3[i]<-1

### we explore the column

              for (l in 1:r)
                {
                  if (l!=i  &  m2[l,j]==1)

                    ## we explore the new row

                    {if (m3[l]==0)  {p <- p-1

                                     m3[l]<-1}

                     for (k in 1:s)

                       {
                         if (k!=j  &  m2[l,k]==1)
                           {
                                        # if (m3[l]==0)  {p <- p-1}
                             m3[l] <- 1
                             m2[i,k] <- 1
                           }
                       }

                     ## we explore the initial row

                     for (k in 1:s)

                       {
                         if (k!=j  &  m2[i,k]==1)
                           {
                                        #  if (m3[l]==0)  {p <- p-1}
                             m3[l] <- 1
                             m2[l,k] <- 1
                           }
                       }
                   }
                } 

### we explore the row

              for (k in 1:s)
                {
                  if (k!=j  &  m2[i,k]==1)

                    ## we explore the new column

                    {		for (l in 1:r)

                                  {
                                    if (l!=i  &  m2[l,k]==1)
                                      {
                                        if (m3[l]==0)  {p <- p-1}
                                        m3[l] <- 1
                                        m2[l,j] <- 1
                                      }
                                  }

                                ## we explore the initial column

                                for (l in 1:r)

                                  {
                                    if (l!=i  &  m2[l,j]==1)
                                      {
                                        if (m3[l]==0)  {p <- p-1}  		 
                                        m3[l] <- 1
                                        m2[l,k] <- 1
                                      }
                                  }
                              }
                }

            }
        }
    }


  return (m2,p)
                                        #return(p)

}

overlappings <- function (data1,data2,indexmethod="intersection") 
{
data1 <- as.matrix(data1)
data2 <- as.matrix(data2)

a <- compare(data1,data2)

s1 <- a$size1
s2 <- a$size2



if (indexmethod == "intersection")
{

 b<- number3(a$int)
 C<- a$int * b$m2
 r <- nrow(C)
 s <- ncol(C)
 p <- b$p

}


else {if (indexmethod == "simpson") 

        { 
          b<- number3(a$s)
          C<- a$int * b$m2
          r <- nrow(C)
          s <- ncol(C)
          p <- b$p
        }

     }

visited <- matrix(0,1,s)
overlapping <- matrix(0,p,3)

#############  compute the overlappings

index <- 1

for (i in 1:r)
 {
   new <- FALSE

   for (j in 1:s)
      {
        if (visited[1,j] == 0) 
           {
	     if (b$m2[i,j]!=0) 
                {
		  new <- TRUE
                  visited[1,j] <- 1
                  overlapping[index,3] <- overlapping[index,3] + sum(C[,j])

		  overlapping[index,1] <-s1 %*% b$m2[,j]
		  overlapping[index,2] <-s2 %*% b$m2[i,]
                }
           }  	
      }
   if (new == TRUE) {index <- index + 1}
 }

#print(overlapping)

###################   compute Jaccard index
R <-matrix(0,p,1)

for (i in 1:p)
  {
    R[i,1] <- R(overlapping[i,1],overlapping[i,2],overlapping[i,3])
  }

#print(c("R scores",R))

#print(c("R_bar",mean(R[,1])))

#print(c("Total R",sum(R[,1])))


return(overlapping,R)

}

colour.bars <- function (y,colours,width = .40) 
{
###### input : a row or a column with p elements from the intersection 
######         matrix, to compute the colour bars

  y<- as.matrix(y)
  p <- max(dim(y))
  size <- sum(y)
  percentages <- round((y*10000)/size)/100
  heights <- c(0, cumsum(y)/size)
  x <- c(0.25,0.25,0.25 + width,0.25 + width)
  plot.new()

  for (i in 1:p) 
    {
      up <- heights[i+1]
      down <- heights[i]
      y.i <- c(down, up, up, down)
      percentage.i <- paste(percentages[i], "%")
      polygon(x, y.i, fg="black",col=colours[i], lwd = 1)
      lines(x = c(0.5,0.53), y = c(up,up), lwd =1, col = colours[i])
    }


}


R <- function (size.A, size.B, intersection) 
{

#### computing the parameter for the termometer

  union <- size.A + size.B - intersection
  R <- intersection / union

  return(R)

}
compare.intersection <- function (datos1,datos2) {
#######################################################################
####
####  function to compute the matrix of intersections between clusters ('int')
####  and the corresponding matrix with Simpson's index ('s')
####
####  the input is a pair of matrices with 'tam' (number of genes) rows and
####  2 columns, the first one for the ID, the second one for the cluster
####  label. The genes are in the same order in both matrices
####
#######################################################################

  datos1<-as.matrix(datos1)
  datos2<-as.matrix(datos2)


########## check the number of clusters in each partitioning  #############
  
  num1<-max(datos1[,2])
  num2<-max(datos2[,2])
  
  tam<-dim(datos1)[1]
  int<-matrix(0,num1,num2)
  
#########  compute the intersections  #################

  for (i in 1:tam) {
    int[as.numeric(datos1[i,2]),as.numeric(datos2[i,2])]<-int[as.numeric(datos1[i,2]),as.numeric(datos2[i,2])]+1
  }

  return(int)
}

compare.simpson <- function (datos1,datos2) {
#######################################################################
####
####  function to compute the matrix of intersections between clusters ('int')
####  and the corresponding matrix with Simpson's index ('s')
####
####  the input is a pair of matrices with 'tam' (number of genes) rows and
####  2 columns, the first one for the ID, the second one for the cluster
####  label. The genes are in the same order in both matrices
####
#######################################################################

  datos1<-as.matrix(datos1)
  datos2<-as.matrix(datos2)

########## check the number of clusters in each partitioning  #############

  num1<-max(datos1[,2])
  num2<-max(datos2[,2])

  tam<-dim(datos1)[1]
  int<-matrix(0,num1,num2)
  
#########  compute the intersections  #################

  for (i in 1:tam) {
    int[as.numeric(datos1[i,2]),as.numeric(datos2[i,2])] <-int[as.numeric(datos1[i,2]),as.numeric(datos2[i,2])]+1
  }

#######  compute the size of clusters  #########

  size1<-matrix(0,1,num1)
  for (i in 1:tam)
    {size1[as.numeric(datos1[i,2])] <-size1[as.numeric(datos1[i,2])]+1
   }
  
  size2 <-matrix(0,1,num2)
  for (i in 1:tam) {
    size2[as.numeric(datos2[i,2])] <-size2[as.numeric(datos2[i,2])]+1
  }


  
########### compute Simpson's index  ############

  s <-int
  for (i in 1:num1) {
    for (j in 1:num2) {
      s[i,j] <-s[i,j]/min(size1[i],size2[j])
    }
  }

  return(s)
}


construct.1 <- function (data, partition) {
#############################################################################
##### INPUT: data -> clustering to be repartitioned
#####
#####        partitions -> 0-1 matrix showing the new partition
#####
#############################################################################

  dat<-as.matrix(data)
  part<-as.matrix(partition)

 ##############################
 # compute the new partitions #
 ##############################

  r<-dim(part)[1]
  s<-dim(part)[2]

  for (i in 1:r) {
    for (j in 1:s) {
      if (part[i,s-j+1]==1) {
        for (l in 1:s) {
          if (part[i,l]==1 & l<s-j+1) {
            dat <- new.partition(dat,s-j+1,l)
          }
        }
      }
    }
  }

  maximum<-max(dat[,2])
  flag<-matrix(0,maximum,1)

  for (i in 1:maximum) {
    for (l in (1:dim(dat)[1])) {
      if (dat[l,2]==i) {
        flag[i] <- 1
      }
    }
  }

  max2<-maximum
  for (i in 1:dim(flag)[1]) {
    if (flag[i]==0 & max2>i) {
      dat <- new.partition(dat,max2,i)
      max2 <- max(dat[,2])
    }
  }

###################################################
###  output: a matrix with the new partitioning ###
###################################################

  return(dat)
}


number <- function (matriz) {
###################################################
#  function to determine a 0-1 matrix to show     #
#  how initial clusters are mapped                #
#  input: any of the matrices of indices obtained #
#  by the functions compara.int or compara.s      #
###################################################
  
################################################
######## MATRIX OF 1'S POINTING OUT THE MAX'S
################################################

  matriz<-as.matrix(matriz)
  r<-dim(matriz)[1]
  s<-dim(matriz)[2]

  m2<-matrix(0,r,s)

  for (i in 1:r) {
    maxim <-max(matriz[i,])
    if (maxim >0) {
      for (j in 1:s) {
        if (matriz[i,j]==maxim) {m2[i,j]<-1}
      }
    }
  }

  for (j in 1:s) {
    maxim<-max(matriz[,j])
    if (maxim>0) {
      for (i in 1:r) {
        if (matriz[i,j]==maxim) {m2[i,j]<-1}
      }
    }
  }


#############################################
######### DETERMINING p #####################
#############################################

  m3<-matrix(0,r,1)
  p<-r

  for (i in 1:r) {
    for (j in 1:s) {
      if (m2[i,j]==1) {
        m3[i]<-1
        ### we explore the column
        for (l in 1:r) {
          if (l!=i  &  m2[l,j]==1) {
            ## we explore the new row
            for (k in 1:s) {
              if (k!=j  &  m2[l,k]==1) {
                if (m3[l]==0)  {p <- p-1}
                m3[l] <- 1
                m2[i,k] <- 1
              }
            }

	    ## we explore the initial row
            for (k in 1:s) {
              if (k!=j  &  m2[i,k]==1) {
                if (m3[l]==0)  {p <- p-1}
                m3[l] <- 1
                m2[l,k] <- 1
              }
            }
          }
        } 

        ### we explore the row
        for (k in 1:s) {
          if (k!=j  &  m2[i,k]==1) {
            ## we explore the new column
            for (l in 1:r) {
              if (l!=i  &  m2[l,k]==1) {
                if (m3[l]==0)  {p <- p-1}
                m3[l] <- 1
                m2[l,j] <- 1
              }
            }
            
	    ## we explore the initial column
            for (l in 1:r) {
              if (l!=i  &  m2[l,j]==1) {
                if (m3[l]==0)  {p <- p-1}  		 
                m3[l] <- 1
                m2[l,k] <- 1
              }
            }
          }
	}
      }
    }
  }

return (m2,p)
#return(p)
}


new.partition <- function (mat,old,new) {
##############################################################
## function to reasign the labels of clusters               ##
## input:  'mat' is a data matrix (genes + cluster labels)  ##
##         'old' is the old label		            ##
##         'new' is the new label                           ##       
##############################################################

  mat<-as.matrix(mat)

  long<-dim(mat)[1]

  for (i in 1:long) {
    if (mat[i,2]==old) {
      mat[i,2]<-new
    }
  }

  return(mat)
}
